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Executive Summary 

An effective interface element technology has been developed for connecting and 
simulating crack growth between independently modeled finite element subdomains (e.g., 
composite plies). This method has been developed using penalty constraints and allows coupling 
of finite element models whose nodes do not necessarily coincide along their common interface. 
Additionally, the present formulation leads to a computational approach that is very efficient and 
completely compatible with existing commercial software. The present interface element has 
been implemented in the commercial finite element code ABAQUS as a User Element 
Subroutine (UEL), making it easy to test die approach for a wide range of problems. The 
interface element technology has been formulated to simulate delamination growth in composite 
laminates. Th anks to its special features, the interface element approach makes it possible to 
release portions of die interface surface whose length is smaller than that of the finite elements. 
In addition, the penalty parameter can vary within the interface element, allowing the damage 
model to be applied to a desired fraction of the interface between the two meshes. Results for 
double cantilever beam DCB, end-loaded split (ELS) and fixed-ratio mixed mode (FRMM) 
specimens are presented. These results are compared to measured data to assess the ability of the 
present damage model to simulate crack growth. 



1. INTRODUCTION 


With model sharing and large scale analysis activities on the rise, there is an increasing 
need to perform a unified analysis of a structural assembly using sub-structural models created 
independently. These sub-structural models are frequently created by different engineers using 
different software and in different geographical locations, with little or no communication 
between die teams of engineers creating the models. As a result, these models are likely to be 
incompatible at their interfaces, making it very difficult to combine them for a unified analysis of 
the entire assembly. Finite element interface technology has been developed to facilitate the 
joining of independently modeled substructures. 


Unconventional approaches have been employed to connect special elements based on 
analytical solutions to finite element models (Aminpour et al., 1991; Jinping et al., 1991). In 
order to take advantage of parallel computing, Farhat et al. (1991,1992) developed a domain 
decomposition approach. In another work (Maday et al., 1988) non-conforming “mortar” 
elements are employed to connect incompatible subdomains. The finite element interface 
technology developed at NASA LaRC (Ransom et al., 1993, 1997; Housner et al., 1995, 
Aminpour et al., 1995, 1997) and elsewhere (Aminpour et al., 1998) allows the connection of 
independently modeled substructures with incompatible discretization along the common 
boundary. This approach has matured to a point that it is now very effective. However, because 
this form of the interface technology utilizes Lagrange multipliers to enforce the interface 
constraint conditions, the resulting system of equations is not positive-definite. Recently, an 
alternative formulation for the finite element interface technology based on Lagrange multipliers 
has been developed (Aminpour et al., 2000). The alternative approach recasts the interface 
element constraint equations in the form of multi-point constraints. This change allows an easier 
implementation of the formulation in a standard finite element code and alleviates the issues 
related to the resulting non-positive definite system of equations. The method seems to provide 
reliable results, but the formulation of the interface method is still quite complicated. 


A possible remedy for these shortcomings is to modify the existing hybrid variational 
formulation of the interface element by enforcing the interface constraints via a penalty method 
as opposed to the current Lagrange multiplier approach. The primary consequences of this 
modification will be (i) a simple formulation that is easily implemented in commercial finite 
elements codes, (ii) a positive-definite and banded stiffness matrix and (iii) a reduced number of 
DOFs, since the Lagrange multiplier degrees of freedom (DOFs) will not be present. Thus, the 
penalty approach should greatly enhance the computational efficiency of the interface element 
technology. 


From an accuracy point of view, the penalty method enforces the constraints only 
approximately, depending on the value of the penalty parameter chosen, while the Lagrange 
multiplier approach enforces the constraints exactly. The penalty method interface approach was 
recently attempted using a single global value of the penalty parameter to enforce all constraints 
(Cho et al., 1998). This study demonstrated the validity and the effectiveness of the penalty 
approach in an interface element However, there is need for specific guidelines regarding the 
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selection of an appropriate value of the penalty parameter, especially when the substructures to 
be connected have different material and/or section stiffnesses. A criterion for choosing the 
penalty parameter in the framework of the interface element under investigation has been 
developed by the authors (Pantano and Averill, 2002a-2002b). 

Though the penalty-based interface element was originally developed to “connect* ' two 
regions of a finite element mesh, it can also be used effectively to simulate crack growth, such as 
delamination. Initiation and evolution of the delamination may be predicted by a fracture 
mechanics approach or by interlaminar strength, introducing an interface constitutive law 
between the layers. 

The fracture mechanics approach uses strain energy release rate G as a parameter for 
assessing delamination, initiation and growth. In the fracture mechanics approach, strain energy 
release rate per unit area delaminated is evaluated and compared with the critical strain energy 
release rate G c . Thus, much attention has been given to correctly evaluating strain energy release 
rate. Pradhan and Tay (1998), Kaczmarek et al. (1998), Rinderknecht and Kroplin (1997), Wang 
and Raju (1996), Raju et al. (1996), Hwu et al. (1995), Hitchings et al. (1996), developed finite 
element approaches for the analysis of delamination growth in composite plates based on the 
virtual crack closure technique (VCCT) to compute the total energy release rate G, including 
contribution by all the three modes. 

In the strength of materials approach, the local state of stress at the interface is compared 
with relevant interface strengths. Use of finite element analysis makes the strength of materials 
approach attractive, as the stresses can be evaluated quickly and efficiently, and interl aminar 
stresses can be easily compared with the measured strengths. However, unless delamination 
initiation is the only concern, failure criteria usually combine strength of materials features with 
fracture mechanics ones. Since the cohesive zone can still transfer load after the onset of 
damage, a softening model is required that describes how die stiffness is gradually reduced to 
zero after the interfacial stress exceeds the interlaminar tensile strength. Reliable prediction of 
the softening behavior can be obtained by relating the work of separation to the critical value of 
the strain energy release rate G c . At a given point on the interface, the area under the stress- 
relative displacement curve is equal to G c . This approach to the determination of the 
delamination growth has been adopted by Reedy et al. (1997), Davila et al. (2001), Mi et al. 
(1998), Chen et al. (1999), Alfano et al. (2001), Lammerant and Verpoest (1996), Schellekens 
and de Borst (1993), Schipperen and Lingen (1999), and Moorthy and Reddy (1999). Interface 
elements are introduced to connect the individual plies of a composite laminate, but the way this 
connection is realized can differ. Two approaches can be identified: point to point interface 
elements acting like springs which connect pairs of nodes, and continuous interface elements that 
connect pairs of two or three dimensional finite elements. 

In the work of Allix and Corigliano (1999), Point and Sacco (1996), Ladeveze (1992), 
Bottega (1983), a constitutive law for the interface material, able to handle the delamination 
phenomenon, has been obtained starting from an adhesion model and is based on the definition 
of a damaged strain energy density in the interface layer. Chakraborty and Pradhan (1999) 
performed a fully 3D finite element analysis to study delamination at the interface of 
graphite/epoxy and glass/epoxy laminates with broken central plies. Their methodology for 
predicting delamination employs simultaneously strength and fracture mechanics approaches. 
Joo and Sim (1994), Ko et al. (1992), Mohammadi et al. (1998), Zhao et al. (1999) and Chang 
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and Springer (1986) predicted delamination initiation by a strength approach very similar to the 
one of Chakraborty and Pradhan (1999). A three-dimensional finite-deformation cohesive 
element and a class of irreversible cohesive laws was developed by Ortiz and Pandolfi (1999) 
and Pandolfi et al. (1999). Moes et al. (1999) presented an innovative finite element able to 
account for crack growth inside the element. The approach adopted for delamination initiation 
and growth is based upon the maximum circumferential stress criterion. Stresses are computed 
by the classical fracture mechanics equation for the stress distribution at the crack tip. This finite 
element type of approach, also referred to as singularity element approach, was also used by 
Aminpour and Hosapple (1991) and Lee and Gao (1995). 

This report will briefly review the penalty-based interface element technology and 
subsequently introduce a new application of the interface element for predicting delamination 
crack growth in laminated structures. 


2. GENERAL DESCRIPTION OF THE INTERFACE ELEMENT 


Consider two independently modeled subdomains Q, and Q 2 as shown in Figure 1(a) 
and 1(b), respectively, for a 2D and for a 3D geometry. The two substructures are connected to 
each other using an interface element acting like “glue” at the common interface. The interface 
element is discretized with a set of nodes that are independent of the nodes at the interface in 
subdomains Q, and 0 2 . The coupling terms associated to the interface element are arranged in 

the form of a “stiffness” matrix and assembled with the other finite element stiffness matrices as 
usual. 


The nodal displacements of die sub-domain Q, are identified by q° and q\ . The 
superscript o marks the degrees of freedom (DOFs) that are not on the interfaces, while i denotes 
DOFs that are on the interfaces. The interface displacement field Uj of the sub-domain Q ,• is 
expressed in terms of the unknown nodal displacements q\. The displacement field V is 

approximated on the entire interface element surface in terms of unknown nodal displacements 

(h- 




V = Tq, 


( 1 ) 


where Nj can be the matrices of linear Lagrange interpolation functions and T is a matrix of cubic 
spline interpolation functions. In the penalty interface method the displacement continuity 
constraint is imposed in a least squares sense using two vectors of penalty parameters y and 75 . 
Thus the total potential energy of the system assumes the form: 


* = *o, + *h, +-7 \( V ~ u i) 2ds+ ^r j( V u if * 


( 2 ) 


The equilibrium configuration is found by taking the first variation of it respect to all the 
DOFs, but not the vectors of penalty parameters y and 75 , which are predetermined constants. 


Sn\ 




= 0 


(3) 


4 


The global system of equations of the penalty hybrid interface method assumes the 
following form: 


(4) 


where: 
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This is a symmetric, banded and positive definite (after boundary conditions are imposed) 
global stiffness matrix. The “stiffness” matrix and generalized vector of unknown displacements 
associated with the interface element can be defined as: 


( 6 ) 


For a detailed description of the interface element formulation, see Pantano and Averill 
(2002a-2002b). 
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3. DETERMINATION OF THE PENALTY PARAMETERS 


In the penalty method, the displacement continuity constraint is imposed through penalty 
parameters, a set of predetermined constants. The FE solution obtained using this method is 
approximate, with its accuracy depending on the value of the adopted penalty parameters. It is 
known that the penalty parameter should depend on the material and/or geometric properties of 
the two sub-regions being joined. Further, there is a relationship between the penalty parameter 
and the Lagrange multiplier that enforces a given constraint. The Lagrange multiplier method 
imposes the continuity constraint exactly; thus it defines the upper limit to the accuracy of the 
penalty method. Knowledge of the correct solution for simple model problem facilitates relating 
the value of the penalty parameter to the geometrical and material properties of the model under 
consideration. In our pursuit of the proper penalty parameter values, a variety of one and two- 
dimensional problems have been studied with both the Lagrange multiplier method and the 
penalty method. The types of finite elements that have been investigated are: conventionally 
formulated and reduced integrated Timoshenko beam elements, plane stress quadrilateral 
elements and plate elements based on the first order shear deformation theory (FSDT), or 
Mindlin plate theory. For each finite element formulation, different penalty parameters are 
associated to the various nodal DOFs. For example, the Timoshenko beam element has three 
independent nodal DOFs: the axial displacement u, the transverse displacement w and the 
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rotation y/. Thus, three different penalty parameters y M , y w and y ¥ are employed to enforce the 
interface continuity constraints on the DOFs u, w and y/. An independent choice of the penalty 
parameters is necessary since each degree of freedom can be related differently to the material 
and geometric properties of die finite element model. 


The methodology adopted in finding the relations will now be described. First the most 
common load cases for the FE type under consideration are applied separately to a simple model 
of one or two elements. The formulations and solutions are obtained using both the Lagrange 
multiplier method and die penalty method. The displacement solutions from the two methods are 
compared individually for each degree of freedom. The ratio between the two solutions is 
expressed in the form: 


^penalty 
u Lagrange 



(7) 


where /=/ (element geometric properties, material properties, and loads). 

Once this simple expression has been identified, the penalty parameter y is set equal to: 
y-Pf . Then, the ratio between the solutions becomes independent of material and geometrical 
properties: 


penalty 1 

—j = 1 + — 

Ugrange p 


( 8 ) 


The accuracy of the solution depends directiy on the value assigned to the parameter p. 
The degree of precision of the solution cannot be indefinitely increased, since round off 
amplification error would rise. However, once a reasonable compromise between constraint 
representation error and the round off error has been evaluated, a value of p can be identified that 
is able to produce the same level of accuracy for every combination of material and geometrical 
Drooerties. 

A A 

For a detailed description of the determination of the penalty parameter values, see 
Pantano and Averill (2002a-2002b). 


4. INTERFACE TECHNOLOGY FOR MODELING DELAMINATION 


The interface element technology (Pantano and Averill, 2002a-2002b) was originally 
developed for the purpose of connecting independently modeled subdomains. By taking the 
converse view, it should be possible to use this technology to disconnect two regions of a model 
during crack growth. In this context, the present interface element approach would have several 
advantages over the conventional FE one. For example, the present approach is able to release 
arbitrary portions of the interface that are smaller than the finite element length, thereby 
determining the crack advance. This can be achieved by changing the extreme values of the 
interval of integration of the interface element or by reducing the value of the penalty parameter 
for a part of that interval. Within each interface sub-region it is possible to evaluate forces at the 
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interface and to reduce the value of the penalty parameter as needed. Thus, it is possible to 
overcome a limitation common to the delamination techniques found in the literature, which 
require delamination growth to be simulated in a discretized form by releasing nodes or elements 
oftheFE model. 

The damage model implemented in die developed interface element is one that is 
commonly adopted (Reedy et al., 1997; Davila et al., 2001; Mi et al., 1998; Chen et al., 1999; 
Alfano et al., 2001; Lammerant and Verpoest, 1996; Schellekens and de Borst, 1993; Schipperen 
and Lingen, 1999; Moorthy and Reddy, 1999). It mixes features of strength of materials 
approaches and fracture mechanics. A bilinear softening model has been implemented in the 
present model (see Figure 2). In single-mode delamination, as the load is progressively 
increased, the relative displacement S the two subdomains FE meshes grows proportionally 
according to the value of the penalty stiffness y. When So is reached the stress is equal to the 
interfacial tensile strength a}, the maximum stress level possible. For higher relative 
displacements the interface accumulates damage and its ability to sustain stress decreases 
progressively. Once S exceeds Sf the interface is fully debonded and it is no longer able to 
support any stress. If die load were removed after So has been exceeded but before Sf has been 
reached, the model would unload to the origin. For example, if after reaching point K (see Figure 
2) the load is reduced, the model unloads along the line KO. If the load is reapplied, the stress 
grows with the relative displacement along the same line KO. 

This behavior is obtained by an effective reduction in the penalty stiffness y. In the 
present model, a new parameter D is introduced in order to signify the damage accumulated at 
the interface: 


o=(\-D)y8 (9) 

Thus D is a damage parameter, whose initial value is zero. D starts growing when 8 > So 
and reaches the value 1 when S > Sf. 

The vaiue of D is computed from geometry to be: 


D(S) = 


S(S F -S B ) 


( 10 ) 


The interfacial constitutive model is entirely defined when any two of the following 
properties are known: G c , o t . So and Sf, where G c is the critical strain energy release rate, which 
is equal to the area under o-8 curve in Figure 2. The following two relations exist among these 
parameters: 

( 11 ) 

4 =— ( 12 ) 

r 

The bilinear interface model is applied to a sub-region of the interface between the two 
meshes; the smaller is the length of each sub-region the higher is the accuracy of the prediction. 
A conventional implementation of the discussed damage technique requires the model to be 
applied along the length of one finite element, wherein the crack can advance in a discrete way 
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only by failing one element at a time. Both limitations necessitate the use of a refined finite 
element mesh. Using the interface element previously presented, the damage evolution scheme 
in our model is effectively mesh-independent, wherein continuity between sub-regions of the 
interface that are much smaller than the finite element length can be released. 

Note that different softening interface laws can be implemented in the present approach. 
In Figure 3 some of the most commonly used damage models are reported. If a damage model 
different from die bilinear one is adopted, the expression (10), which allows die evaluation of 
damage parameter D as a function of the relative displacement S, must be modified. The changes 
in the relation (10) will easily follow from the new geometry of the softening interface model. 

Thus, the present interface model can be applied to a desired fraction of the interface 
between the two meshes. If we divide the interface element into a given number of intervals w, 
each of them will obey the rules of the failure model independently from the others. This 
corresponds to changing the total potential energy of the system in the following way: 

* = *-n,+*n 2 +^Z(l-A)r,£(^-« 1 ) 2 * + Z( 1 -A)r 2 £ j (^-«0 2 * (13) 


where Z), is the damage parameter associated with the interval i, and the interval i is defined over 
the range ( Z,_y, Z, ). Z, is the value of the interface coordinate Z at the end of the i* interval. The 
value of the relative displacement 5 is evaluated at the center of the interval i. By allowing crack 
advance in a more continuous way, a higher accuracy of the simulation can be achieved. 


To illustrate this concept, consider the two incompatible meshes shown in Figure 4. They 
are joined by two interface elements whose length equals that of five conventional elements of 
the top mesh and four elements of the bottom mesh. Two forces F applied at the tip are 
responsible for effecting a mode I stress field at the interface. The interface element at the crack 
tip is shown in Figure 4(a) as divided in ten intervals. The intervals do not need to match any of 
the nodes in the upper or lower mesh, but in this example some of them coincide. Moreover, the 


number of intervals in the interface element is a parameter that can be changed as desii 


Following the progression in Figure 4, a simulation of the delamination growth is achieved by 
releasing portions (intervals) of the interface element. In Figure 4(b), the first interval is failed. 
The portion of the interface element not released still applies its constraint to the lower and 
upper elements next to the crack tip. 


In Figure 4(c) the second interval is failed, which determines the complete release of the 
element on the upper mesh near the crack tip. The element on the lower mesh moves downward 
too, but being still held in part, the movement is small. The next advance in the crack length, 
Figure 4(d), frees die lower element and only partially the upper. Figure 4(e) shows the effect of 
releasing yet another interval. Similarly, when all the intervals are failed die FE model behaves 
as though the first interface element is not present. Then, the next interface element starts failing. 
In this example, intervals whose length is half of the smaller finite element extension have been 
used. Dividing the interface element into more intervals or reducing its length would improve the 
accuracy of the model. 


By progressively reducing the dimension of the intervals in which the interface element 
is divided, the results will converge to a values that will not change noticeably as the intervals 
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length are further reduced. Usually, this kind of convergence study do not require several 
simulations since the rate of convergence is typically high. Moreover, experience in choosing the 
intervals length developed by working with similar simulations can be easily applied to new 
calculations. It needs to be underlined that intervals dimensions smaller than the minimum 
element size along the crack should not be used as the accuracy will quickly be reduced. In the 
numerical results convergence of the solution with the number of intervals is investigated for the 
double cantilever beam DCB specimen. 

The other essential convergence study, which is commonly necessitated in FE approaches 
to the simulations of crack growth, is based on investigating the model behavior as a function of 
the number of increments in which the given load/displacement is progressively applied. At the 
beginning of each increments, or load step, the relative displacement 8 at the center of the 
interval is obtained from the FE calculation for all the intervals in which die interface element is 
divided. If the increment in the applied load or displacement is too high, the bilinear softening 
model cannot work correctly and die results will be incorrect For example, if at the beginning of 
the increment the relative displacement 8 for the interval at the crack tip is slighdy smaller than 
So while at the end of the increment it is much higher, for the entire load step the damage 
parameter D maintain its zero value and the stress at the interface will considerably exceed the 
interfacial tensile strength a,. The softening will be applied only in the next increment, with the 
additional flawed consequence that when the interval will fail the area under 0-8 curve will be 
higher than the critical strain energy release rate G c . Thus, as for the great majority of the FE 
approaches to the simulations of crack growth, a convergence study on the number of increments 
should always be performed. In the numerical results convergence of the solution with the 
number of increments is reported for the double cantilever beam DCB specimen. 

For a more detailed description of the interface element technology for simulating crack 
growth between independently modeled finite element subdomains, see Pantano and Averill 
(2002b). 

5. Mixed Mode Failure 


In defining the softening model, the following properties are needed: the interlaminar 
tensile and shear strengths T and S, the penalty parameter y, and the critical strain energy release 
rates Gj c and Gu c . For simplicity, we assume herein that material behavior is the same for both 
tensile and compressive loading. 

At a given load increment the finite element solution allows us to compute 8 and 8 for an 
interval of the interface element. We also know: 

< 5,0 = — 04 ) 
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( 15 ) 
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cr z = r s z (16) 
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( 17 ) 


r « = r 8 

A quadratic interface failure criterion can be written in the following forms: 
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If the condition in (19) is not satisfied, no action needs to be taken. Otherwise, we assume the 
following ratio: 


'a' 1 

^-=c 




( 20 ) 


to be the same as it was when the failure condition (19) was satisfied first (see Figure 5). 

If die load steps are small, then the assumption should be valid. Then, from geometry, it 
is possible to determine the value of the relative displacements 4’ and 8 corresponding to point 
Fin Figure 5. 
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The interfacial constitutive models are updated for both delamination modes I and II by setting 
8 x 0 = 8 x’ and 8 x 0 = 8 ' (see Figures 6 and 7). 

The interlayer tensile strengths T and S are updated accordingly: 

r=y 2 S' z o (23) 

S'=r x - 8 ' x o (24) 

Note that the following inequalities hold. 


‘V^O. S Z ,'<S M , 8 x ^ 8 x0 ' , 8 X > 8 X ,' 


(25) 
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The quadratic interaction criterion predicts final failure to be reached when the following 
condition is satisfied: 



In a similar manne r as before, we assume the ratio between (Gj/GjjtJ and (G/Giq) does not vary 
as the work of separation grows, as shown in Figure 8: 



(27) 


Research on the specimens commonly used for delamination studies shows that, for a given 
configuration (geometry and loads), the ratio between the strain energy" release rates for modes I 
and II, Gj/Gn, does not change appreciably during the entire test (e.g. Choi et al., 1999). This fact 
provides a valid foundation to our assumption. Thus, it is possible to determine the value of the 
strain energy release rate G{ and Gu corresponding to point F in Figure 8. 

From geometry it is possible to determine the value of the (G/Gn) ' at F. 



(28) 


In order to evaluate this expression we need to determine (Gj/Giq), (actual value of Gj divided by 
Gic), see Figure 9. 
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Similarly, we have: 
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Now, the updated state of the interfacial constitutive models can be completed for both modes I 
and II by setting Gj c ' — Gf and Gnc = Gjj. 

Figures 9 and 10 show the final form of the interfacial constitutive models. As can be noticed, 
the models have different penalty and damage parameters. In this a way, the maximum freedom 
is allowed for modeling delamination modes I and n. 


6. Friction Model 

6. 1 Evaluation of the force at the interface 


It is known that the Lagrange multiplier X is related to the penalty parameter yby: 

F = A = y* J(v-u)ds (33) 

s 

Where F is die force per unit length required to hold the two regions together. It is important to 
notice that we are not restricted to compute die force on the entire interface, but rather it can be 
evaluated for a portion of the interface length. This can be performed by changing the extreme 
values of the interval of integration. 

Finally, computation of the interface force does not depend on the compatibility of the 
interface meshes. Even if die discretizations of the adjacent layers, or layer and skin-stiffener, 
are different we are still able to choose at will the interface length along which to evaluate 
forces. 


6.2 Implementation of Friction in the Interface Element 


The friction model can be applied to interface elements after complete failure or to 
interface elements whose only purpose is to avoid overlapping and to enforce friction. 

For each of the intervals in the interface, the normal force F„ at the interface can be 
calculated in terms of the normal relative displacement 8„ : 

< 34 > 

Assuming Coulomb friction, the tangential force F, that the interface element needs to 
generate to simulate the friction phenomenon is: 
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(35) 


F,=\r,('-D,)\6,<Is = TF. 


where ris the friction coefficient and 5, is the tangential relative displacement. The parameter D t , 
which before failure was used as a damage parameter, is now employed as a scale factor able to 
reduce the value of the penalty parameter for tile tangential DOF. Assuming the value of D, that 
makes the equality (35) to hold true, the right amount of friction will be generated by the 
interface element 


The required damage parameter D] related to the tangential relative displacement 5, can 
be determined from the following relation: 



s 


(36) 


6. 3 Numerical Test of the Friction Model 


A simple test has been performed to verify the reliability and the accuracy of die friction 
model implemented into the interface element 

The loading, the boundary conditions and the geometry for the test problem are 
illustrated in Figure 11. Two finite element meshes are connected along the common boundary 
using one interface element. All the finite elements are rigid and their thickness is 0.5m. 

During a quasi-static analysis of 400 load increments, a 1mm tangential displacement is 
gradually imposed to the nodes of the upper mesh, as shown. At the same time, a pressure load 
of 1 MPa is applied to the top the upper mesh and is kept constant through the entire analysis. 
According to the Coulomb law F, — t-N, since a coefficient of friction equal to r = 0.1 is 
assumed and the contact area is 1 m 2 , as expected total reaction force at the interface nodes is: 

F t = tF„ = to A = 0.1 • 1 ,e6^r ■ 1 jn 2 = 1 ,e5N (37) 

m 

In the Figure 12 is reported the reaction forces versus the tangential displacement 
produced by the FE analysis for the interface nodes 1, 2 and 3. For each of the 400 increments in 
die displacement value, the reaction force at the central node is recorded to be a constant value of 
0.5e5N while on the other two nodes a constant value of 0.5e5N is obtained. These sum to the 
analytical value of 1 .e5 N. 
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7. NUMERICAL RESULTS 


7.1 Double Cantilever Beam Test #7 


The loading, the boundary conditions and the geometry for the double cantilever beam 
DCB specimen are illustrated in Figure 13. The DCB test is recognized as pure mode I test. The 
properties assumed for the beam material are: 2s//=130 GPa, E 22 - E 33 - 8 GPa, G 12 = 6 GPa, 
v- 0.27. The properties of the DCB specimen interface are: G/ c = 257 N/m and <x, = 20 MPa. 

Three finite element models of varying mesh refinement (120x2, 300x8, 600x8) have 
been used to verify the capabilities of the present approach. In Figure 14, original and deformed 
models of the DCB test specimen are shown for the coarsest mesh. Two independent meshes 
compose the finite element models of the upper and lower part of the beam; they are joined by 
several interface elements. Each interface element connects four finite elements: two on the 
bottom mesh and two on the top. Experimental results used to validate the present delamination 
approach have been reported by Chen et al. (1999). A plot of the reaction force as a function of 
the applied end displacement is shown in Figure 15. Results from the experiment are compared 
to analytical ones for all the finite element models. 

Note that the total number of elements in the meshes is indicated, for example two 120x1 
meshes joined by interface elements compose the 120x2 mesh. The predictions from the 120x2 
mesh contain many local “bumps”. Mi et al. (1999) have analyzed this phenomenon, concluding 
that coarse meshes can induce these “false instabilities”. As the mesh is refined, the predicted 
response agrees very well with that measured experimentally. 

As discussed previously, the damage technique implemented in our model allows 
portions of the interface, intervals, much smaller than the finite element length, to be released. In 
Figure 16 convergence of the solution with the number of intervals is investigated for the 300x8 
mesh. It can be noticed that as the number of intervals increases the accuracy of the results 
improves. Figure 17 illustrates the model behavior as a function of the number of increments for 
the 300x8 mesh. Convergence of the solution to the experimental one is achieved by increasing 
number of intervals. 


7.2 Double Cantilever Beam Test #2 


To further verify the ability of the interface to accurately simulate delamination growth in 
composite materials, a second double cantilever beam DCB test was performed. Material data 
and experimental results are taken from a different source (Alfano et al., 2001) with respect to 
the first DCB test 

The loading, the boundary conditions and the geometry for the double cantilever beam 
DCB specimen are illustrated in Figure 18. The properties assumed for the beam material are: 
En = 126 GPa, E22=E 3 3=7.5 GPa, G/r=4.981 GPa, v =0.281. The properties of the DCB specimen 
interface are: G Ic = 263 N/m and a, — 57 MPa. 
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A plot of the reaction force as a function of the applied end displacement is shown in 
Figure 19. The finite element model has a mesh of 300x8 elements and the displacement 
boundary conditions are applied in 400 increments. The predicted results compare well with the 
experimental results. 


7.3 End-Loaded Split (ELS) Beam Test 


To characterize the mode II delamination, shearing mode, the end-loaded split (ELS) 
specimen is often used. The loading, the boundary conditions and the geometry for the end- 
loaded split (ELS) specimen are illustrated in Figure 20. 

The properties assumed for the beam material are: £//=130 GPa, En~ E 33 = 8 GPa, Gn = 
6 GPa, v= 0.27. The properties of the ELS specimen interface are: Gjj c = 856 N/m and cr, = 48 
MPa. As in the DCB test specimen models, two meshes compose the finite element models of 
the ELS specimen; they are joined by several interface elements. For the results shown, a 300x8 
finite element mesh has been utilized. 

In Figure 21, original and deformed models of the ELS test specimen are shown. 
Experimental results used to validate the present delamination approach have been reported by 
Chen et al. (1999). A plot of the reaction force as function of the applied end displacement is 
shown in Figure 22. Convergence of the solution as the number of increments increases is 
demonstrated. The comparison with experimental results is again very good. 


7.4 Fixed-Ratio Mixed Mode (FRMM) Test 

To test the mixed mode I+n delamination capabilities of the interface model, the Fixed- 
Ratio Mixed Mode (FRMM) specimen has been considered. The loading, the boundary 
conditions and the geometry for the FRMM specimen are illustrated in Figure 23. 

The properties assumed for the beam material are: £^=130 GPa, E 22 — E 33 = 8 GPa, Gn = 
6 GPa, v= 0.27. The properties of the ELS specimen interface are: G/ c = 257N/m, Gu c - 856N/m 
and o, = 48 MPa. As in the DCB and ELS test specimens models, two meshes compose the finite 
element models of die ELS specimen; they are joined by several interface elements. For the 
results shown, a 300x8 finite element mesh has been utilized. 

In Figure 24, original and deformed models of the ELS test specimen are shown. 
Experimental results used to validate die present delamination approach have been reported by 
Chen et al. (1999). A plot of the reaction force as function of the applied end displacement is 
shown in Figure 25. The predictions agree well with the experimental results, demonstrating that 
the novel mixed mode damage model implemented in the interface element is able to correctly 
predict delamination growth in a mixed mode FE simulation. 
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8. CONCLUSIONS 


A new interface element technology has been presented for predicting crack growth in 
laminated structures. This interface method is capable of joining and simulating crack growth 
between independently modeled finite element subdomains (e.g., composite plies). The interface 
’ element approach makes it possible to release sub-regions of the interface surface whose length 
is smaller than that of the finite elements, thereby allowing for a mesh-independent tracking of 
the crack front A novel approach for modeling the crack growth in mixed mode I+II conditions 
has been developed. Results for double cantilever beam DCB, end-loaded split (ELS) and fixed- 
ratio mixed mode (FRMM) specimens indicate that the method is capable of accurately 
predicting delamination growth. 
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Figure 3. Interfacial constitutive damage models commonly adopted, 
where pp refers to , pro refers to ... , tin refers to the bilinear 
model used herein, and reg refers to ... . 
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Figure 11. Loading, boundary conditions and geometry for the friction 
test problem, a) Exploded view, b) assembled view. 
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Figure 12 . Test Friction - Graph Force-Displacements. 
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Figure 13. Geometry and boundary conditions for the DCB test 
specimen. 



Figure 14. Original and deformed (5mm opening) models of the DCB 
test specimen using a 120x2 mesh. 
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Figure 15. Force vs. displacement results of experimental and numerical 
DCB test 
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Figure 16. Force vs. displacement results of DCB test. Convergence of 
the solution with intervals number. Mesh 300x8. 
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Figure 17. Convergence of the solution with the number of increments. 
Mesh 300x8. 
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Figure 19. Force vs. displacement results of experimental and numerical 
DCB test #2 
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Figure 20. Geometry and boundary conditions for the ELS test 
specimen. 



Figure 21. Original and deformed (30mm tip displacement) models of 
the ELS test specimen for a 300x8 mesh. 
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Figure 22. Force vs. displacement results of experimental and analytical 
ELS test for varying number of loading increments. 
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Figure 23. Geometry and boundary conditions for the FRMM test 



Figure 24. Original and deformed (20mm tip displacement in absence of 
delamination) models of the FRMM test specimen. 
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Figure 25. Force vs. displacement results of experimental and numerical 
FRMM test. 
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